*===========================================================================		
// Density calculations
capture program drop calculate_density
program calculate_density
syntax [anything(name=var)], [knn(real -1) levels(name) cont]

	if `knn' == -1 {
		mata density = ((M*proximity):/(J(Nix,Npx,1)*proximity))
	}
	else {
		mata simi = J(Npx,Npx,0)
		forval i = 1/$Np { 			
			mata temp = sort(proximity[.,`i'],-1)
			mata simi[.,`i'] = proximity[.,`i'] :* (proximity[.,`i'] :>= temp[`knn'])
		}
		mata weight = simi:/(J(Npx,Npx,1)* simi)
		
		if "`cont'"~="" {
			noi display "	: Continuous"
			noi di "`levels'"
			if("`levels'" == "") {
				mata density = RCA * weight
			}
			else {				
				mata density = `levels' * weight 
			}
		}	
		else {
			noi display "	: Discrete"
			mata density = M * weight
		}	
	}
end	
*===========================================================================
